function NotchBeamDamage
clc
%close all
addpath('../');
addpath('../../');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% read mesh
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
projDir = '/home/hbui/kratos_janosch';
%coarse (structure) mesh
projName = 'notchbeam.gid';
fn = [projDir '/' projName '/' 'mesh1.msh'];

mg = fem2.gid.GidMeshGenerator(fn);
mesh = mg.getMesh2d();
mesh.setNumberOfFields(2);
%set id for nodes
for i=1:mesh.countNodes()
    mesh.getNode(i-1).setId(i);
end
%set id for faces
for i=1:mesh.countFaces()
    mesh.getFace(i-1).setId(i);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% conditions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% lambda = [linspace(0,1,10) linspace(1,1.5,11)];
du = 2;
%lambda = [linspace(0,1,30)]; %work with du = 2
lambda = [linspace(0,1,60)];
% du = 4;
% lambda = [linspace(0,1,100)];

% du = 2.376;
% lambda = [linspace(0,1,30)];
meshUtil = fem2.MeshUtilities(mesh);
c1 = fem2.Constraint([0 0]);
c2 = fem2.Constraint([0 1]);
c3 = fem2.Constraint([0 1]);
c3.setValue(0,du);
edge1 = meshUtil.seekNode(1,0,0);
edge2 = meshUtil.seekNode(1,0,-2);
mesh.setConstraint(c2,edge1);
meshUtil.seekNode(0,0).setConstraint(c1);
mesh.setConstraint(c3,edge2);
% plotMesh(mesh);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% element declaration
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
M = fem2.Model(mesh);
nf = mesh.countFaces();
E = 2;
nu = 0.3;
ft = 0.91;
e0 = ft/E;
ef = 50*e0;
DamageLaw = 'exponential';
% EquivalentStrainMode = 'euclidean';
EquivalentStrainMode = 'strain norm';
C = C_plane_stress(E,nu);
for i=1:nf
    e = fem2.matlab.T3MatlabElement(mesh.getFace(i-1));
    M.addElement(e);
    %create struct to contain element state variables
    es = createElementContainer();
    es.constitutiveMatrix = C;
    es.elasticConstitutiveMatrix = C;
    es.e0 = e0;
    es.ef = ef;
    es.kappa = e0;
    es.kappa_old = e0;
    es.kappa_old_old = e0;
    es.DamageLaw = DamageLaw;
    es.EquivalentStrainMode = EquivalentStrainMode;
    elements(i) = es;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solve
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:edge1.size()
    edge1nodes(i) = edge1.get(i-1).getId();
end

gid = fem2.gid.GidPost(projDir,projName);

% postProcessInfo.infoList = {'resistant force',...
%     'strain plot for softening'};
% h = figure();
% postProcessInfo.infoParam = {edge1nodes,{h,[7]}};

% postProcessInfo.infoList = {'resistant force',...
%     'write damage at Gauss points to gid'};
% postProcessInfo.infoParam = {edge1nodes,gid};

postProcessInfo.infoList = {'resistant force'};
postProcessInfo.infoParam = {edge1nodes};

[elements,output] = solving_strategy_static_3(elements,M,...
    lambda,50,1e-6,postProcessInfo);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% post process
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
res_node_2(:,:) = -output.resistant_force(:,2,:);
delta = lambda*du;
plot(delta,res_node_2);

rmpath('../');
rmpath('../../');
